
# Code to create Figure 1:  projected changes in temperature and precipitation over US corn growing regions
# Also creates companion figure for Africa (Appendix Fig A1)

# Fig 1 for US
setwd("C:\\Users\\Elisa Cascardi\\Desktop\\BDLMS_Replication\\")  #set working directory to the directory where the replication .zip was unzipped. Include final slash. You may need to replace single backslash with double backslash depending on operating system.

load('data\\additionalForFigures\\clim_mod_chg_US.Rdata')
modelnames=read.csv("data\\additionalForFigures\\modelnames.csv")
models=as.vector(unique(unlist(modelnames),na.rm=T))[1:20]  #unique model names
runs=c()  #this will give you which models are in which runs
for (k in 1:length(models)) {
	i=models[k]
	add=c()
	if (i%in%modelnames[,1]) {add=c(add,1)}
	if (i%in%modelnames[,2]) {add=c(add,2)}
	if (i%in%modelnames[,3]) {add=c(add,3)}
	runs[[k]]=add
	}

# now make the figure
#pdf(file="output\\Figure1.pdf",height=6,width=6)
cols=c("grey","grey20","white")
plot(1,type="l",xlim=c(1,8.2),ylim=c(-49,25),xlab="temperature change (deg C)",ylab="precip change (%)",las=1,axes=F)
# par(mar=c(5,4,2,1))
for (i in 1:length(models)) {
	toplot=runs[[i]]
	x=y=c()
	for (a in toplot) {
		x=c(x,chg[[a]][which(modelnames[,a]==models[i]),1])
		y=c(y,chg[[a]][which(modelnames[,a]==models[i]),2]*100)
		}
	if (length(toplot)==2) {lines(x,y,lwd=1,lty=3)}
	if (length(toplot)==3) {
		lines(x[1:2],y[1:2],lwd=1,lty=3)
		lines(x[c(1,3)],y[c(1,3)],lwd=1,lty=3)
		}
	}
for (i in 1:length(models)) {
	toplot=runs[[i]]
	x=y=c()
	for (a in toplot) {
		x=c(x,chg[[a]][which(modelnames[,a]==models[i]),1])
		y=c(y,chg[[a]][which(modelnames[,a]==models[i]),2]*100)
		}
	if (models[i]=="hadcm3") {
		points(x,y,pch=24,bg=cols[toplot],col="black",cex=1.8) 
		} else {
		points(x,y,pch=21,bg=cols[toplot],col="black",cex=1.8)	}
	}
ath=c(-41,-49,-45)
for (i in 1:3) {
	toplot=chg[[i]][,1]
	lowhi=quantile(toplot,probs=c(0.025,0.975))
	toplot=replace(toplot,toplot<lowhi[1],lowhi[1])
	toplot=replace(toplot,toplot>lowhi[2],lowhi[2])
	boxplot(toplot,at=ath[i],horizontal=TRUE,range=0,add=T,boxwex=4,col=cols[i],axes=F)
	}
atv=c(7.5,8.2,7.85)
for (i in 1:3) {
	toplot=chg[[i]][,2]*100
	lowhi=quantile(toplot,probs=c(0.025,0.975))
	toplot=replace(toplot,toplot<lowhi[1],lowhi[1])
	toplot=replace(toplot,toplot>lowhi[2],lowhi[2])
	boxplot(toplot,at=atv[i],range=0,add=T,boxwex=0.3,col=cols[i],axes=F)
	}
abline(h=-38)
abline(v=7.3)
axis(2,at=seq(-30,20,10),seq(-30,20,10),las=1)
axis(1,at=1:7,1:7)
box()
lines(c(5.1,5.5),c(1,6))
text(6,7.3,"HadCM3 A1b",cex=0.8)
s=c("A1b","A2","B1")
text(0.65,ath,s,cex=0.7,pos=4)
text(atv,22,s,srt=90,cex=0.7)
#dev.off()
dev.copy2eps(file="output\\Figure1.eps",height=6,width=6)

#calculation of model range represented by 2 models chosen randomly, as a % of total model range. based on A1b scenario, temperature.  also calculate whether actual ensemble mean falls between two models
x=chg[[1]][,1]
rg=max(x)-min(x)
mn=mean(x)
rgg=c()
inn=c()
for (i in 1:10000) {
	y=sample(x,2)
	rg1=max(y)-min(y)
	rgg=c(rgg,rg1/rg)
	inn=c(inn,mn<max(y)&mn>min(y))
	}
mean(rgg)
mean(inn)


###################################################################################################
# Figure A1 for Africa
###################################################################################################

load('data\\additionalForFigures\\clim_mod_chg_Afr.Rdata')
modelnames=read.csv("data\\additionalForFigures\\modelnames.csv")
models=as.vector(unique(unlist(modelnames),na.rm=T))[1:20]  #unique model names
runs=c()  #this will give you which models are in which runs
for (k in 1:length(models)) {
	i=models[k]
	add=c()
	if (i%in%modelnames[,1]) {add=c(add,1)}
	if (i%in%modelnames[,2]) {add=c(add,2)}
	if (i%in%modelnames[,3]) {add=c(add,3)}
	runs[[k]]=add
	}

#pdf(file="output\\Figure_A1.pdf",height=6,width=6)
#quartz(width=6,height=6)
cols=c("grey","grey20","white")
plot(1,type="l",xlim=c(1,5.7),ylim=c(-14,15),xlab="temperature change (deg C)",ylab="precip change (%)",las=1,axes=F)
for (i in 1:length(models)) {
	toplot=runs[[i]]
	x=y=c()
	for (a in toplot) {
		x=c(x,afrchg[[a]][which(modelnames[,a]==models[i]),1])
		y=c(y,afrchg[[a]][which(modelnames[,a]==models[i]),2]*100)
		}
	if (length(toplot)==2) {lines(x,y,lwd=1,lty=3)}
	if (length(toplot)==3) {
		lines(x[1:2],y[1:2],lwd=1,lty=3)
		lines(x[c(1,3)],y[c(1,3)],lwd=1,lty=3)
		}
	}
for (i in 1:length(models)) {
	toplot=runs[[i]]
	x=y=c()
	for (a in toplot) {
		x=c(x,afrchg[[a]][which(modelnames[,a]==models[i]),1])
		y=c(y,afrchg[[a]][which(modelnames[,a]==models[i]),2]*100)
		}
	if (models[i]=="hadcm3") {
		points(x,y,pch=24,bg=cols[toplot],col="black",cex=1.6) 
		} else {
		points(x,y,pch=21,bg=cols[toplot],col="black",cex=1.6)	}
	}
ath=c(-11,-14,-12.5)
for (i in 1:3) {
	toplot=afrchg[[i]][,1]
	lowhi=quantile(toplot,probs=c(0.025,0.975))
	toplot=replace(toplot,toplot<lowhi[1],lowhi[1])
	toplot=replace(toplot,toplot>lowhi[2],lowhi[2])
	boxplot(toplot,at=ath[i],horizontal=TRUE,range=0,add=T,boxwex=1.8,col=cols[i],axes=F)
	}
atv=c(5.2,5.7,5.45)
for (i in 1:3) {
	toplot=afrchg[[i]][,2]*100
	lowhi=quantile(toplot,probs=c(0.025,0.975))
	toplot=replace(toplot,toplot<lowhi[1],lowhi[1])
	toplot=replace(toplot,toplot>lowhi[2],lowhi[2])
	boxplot(toplot,at=atv[i],range=0,add=T,boxwex=0.25,col=cols[i],axes=F)
	}
abline(h=-10)
abline(v=5)
axis(2,at=seq(-10,15,5),seq(-10,15,5),las=1)
axis(1,at=1:5,1:5)
box()
s=c("A1b","A2","B1")
text(0.75,ath,s,cex=0.7,pos=4)
text(atv,15,s,srt=90,cex=0.7)
lines(c(4,3.5),c(-7,-0.2))
text(4.3,-7.5,"HadCM3 A1B",cex=0.8)
#dev.off()
dev.copy2eps(file="output\\FigureA1.eps",height=6,width=6)


